Monte Carlo simulation with time step quantification in terms of Langevin dynamics 
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For the description of thermally activated dynamics in systems of classical magnetic moments 
numerical methods are desirable. We consider a simple model for isolated magnetic particles in a 
uniform field with an oblique angle to the easy axis of the particles. For this model, a comparison 
of the Monte Carlo method with Langevin dynamics yields new insight in the interpretation of the 
Monte Carlo process, leading to the implementation of a new algorithm where the Monte Carlo 
step is time-quantified. The numeric results for the characteristic time of the magnetisation reversal 
are in excellent agreement with asymptotic solutions which itself are in agreement with the exact 
numerical results obtained from the Fokker-Planck equation for the Neel-Brown model. 



Studies of spin dynamics in particulate systems are 
currently of significant interest, as model systems for un- 
derstanding the thermodynamics of the reversal process. 
Brown m developed a theoretical formalism for thermally 
activated magnetisation reversal based on the Fokker- 
Planck (FP) equation which led to a high energy bar- 
rier asymptotic formula in the axially symmetric case of 
a particle with easy (uniaxial) anisotropy axis coUinear 
with the applied magnetic field. Since then extensive 
calculations P-ffl have been carried out in which im- 
proved approximations were found for the axially sym- 
metric case. Coffey and co-workers 0-0] also derived 
formulae for the non-axially symmetric case, investigat- 
ing also the different regimes imposed by the damping 
parameter a of the Landau-Lifshitz-Gilbert (LLG) equa- 
tion. This work represents an important basis for the un- 
derstanding of dynamic processes in single-domain par- 
ticles. New experimental techniques which allow for an 
investigation of nanometer-sized, isolated, magnetic par- 
ticles confirmed this theoretical approach to thermal ac- 
tivation [Q. 

Unfortunately, the extension of this work to the im- 
portant case of strongly coupled spin systems such as are 
found in micromagnetic calculations of magnetisation re- 
versal is non-trivial, and realistic calculations in systems 
with many degrees of freedom would appear to be impos- 
sible except by computational approaches. These are cur- 
rently of two types: (i) calculations involving the direct 
simulation of the stochastic (Langevin) equation of the 
problem, in this case the LLG equation supplemented by 
a random force representing the thermal perturbations. 
This is referred to as the Langevin Dynamics (LD) for- 
mahsm [|[||, and (ii) Monte Carlo (MC) simulations Q 
with a continuously variable (Heisenberg like) Hamilto- 
nian |l^Jl3]. The LD approach, although having a firm 
physical basis is limited to timescales of the order of a 
few ns for strongly coupled systems. The MC approach 
is capable of studying longer timescales involving rever- 



sal over large energy barriers, but has the severe problem 
of having no physical time associated with each step, re- 
sulting in unquantified dynamic behavior. 

Physically, the dynamic behavior of interacting spin 
systems is a topic of considerable current interest, much 
of this interest being driven by the need to understand 
spin electronic devices such as MRAM. The possibility of 
truly dynamic models of strongly coupled systems would 
seem to be an important factor in the development of a 
fundamental physical understanding. This requires dy- 
namic studies over the whole time range from ns and 
sub - ns to the so-called 'slow dynamic' behavior arising 
from thermally excited decay of metastable states over 
timescales from 10-lOOs and upwards. It is inconceiv- 
able that the LD technique can be used over the whole 
timescale and therefore a truly time quantified MC tech- 
nique is necessary in order to allow calculations over the 
longer timescales of physical interest. Here we propose a 
technique for the quantification of the MC timestep and 
give a supporting argument developed from the fluctu- 
ation dissipation theorem. This argument results in a 
theoretical expression for the timestep in terms of the 
size of MC move, and also gives the validity criterion 
that the MC timestep is much longer than the precession 
time. Comparison with an analytical formula for relax- 
ation in the intermediate to high damping limit is used 
to verify the theoretically predicted relationship relating 
the timestep to the size of MC move. This represents 
an important first step in the process of deriving a theo- 
retical formalism for time quantified MC calculations of 
strongly interacting spin systems. 

We consider an ensemble of isolated single-domain par- 
ticles where each particle is represented by a magnetic 
moment with energy 



E{S) = -dVS^, - ^isB ■ S, 



(1) 



where S_ — i-t/ fJ-s is the magnetic moment of unit length, 
_B = BxX + BzZ represents a magnetic field under an ar- 
bitrary angle ip to the easy axis of the system, d is the 



uniaxial anisotropy energy density and V the volume of 
the particle. Throughout the article we use the mate- 
rial parameters T^ = 8 x 10^^'^m'^, d — A.2 x lO^J/m , 
magnetic moment /i^ = 1.12 x lO^^^J/T. 

The LLG equation of motion with LD R] is 



{l + a^)fi. 



S X (H{t) +aSx H{t) 



(2) 



where 7 = 1.76 • 10"'^^ (Ts) ^ is the gyromagnetic ra- 
tio, Hit) = C{t) — ^, and Q is the thermal noise with 
(C.(t)) = and {Ut)Q{t')) = 5,j6{t - t')2akBTii,/j. i 
and j denote the cartesian components. 

The equation above is solved numerically using the 
Heun method |Q]. Also, it is possible to obtain analyti- 
cally asymptotic solutions for the escape rate which have 
been extensively compared with the exact numerical so- 
lutions from the corresponding matrix form of the FP 
equation for a wide range of parameters and non-axially 
symmetric potentials Q-g]. 

Both of our simulations, MC as well as LD, start with 
the magnetic moments in z-direction. The magnetic field 
has a negative z-component so that the magnetization 
will reverse after some time. The time that is needed for 
the z-component of the magnetization to change its sign 
averaged over a large number of runs {N = 1000) is the 
characteristic time r which corresponds to the inverse of 
the escape rate following from exact numerical solutions 
of the corresponding FP equation. 

For the MC simulations we use a heat-bath algorithm. 
The trial step of our MC algorithm is a random move- 
ment of the magnetic moment within a cone with a given 
size. In order to achieve this efficiently we construct 
a random vector with constant probability distribution 
within a sphere of radius R. This random vector is added 
to the initial moment and subsequently the resulting vec- 
tor is normalized. 

The size of the cone R of our algorithm influences the 
time scale the method simulates. We investigate the in- 
fluence of R on our MC algorithm by varying R and 
calculating r. As usual in a MC procedure the time 
is measured in Monte Carlo steps (MCS). For our cal- 
culation we use a field of \B_\ = 0.2T and an angle of 
■0 = 27° to the easy axis. The resulting energy barrier is 
AE — 8.2 X 10~^^J, the temperature we chose for Fig. n^ is 
AE/ksT = 3.3. As Fig. demonstrates, it is r ~ R'^. 
This dependence can be understood by considering the 
moments as performing a random walk where R is pro- 
portional to the mean step width. Having understood 
that the MC time can be set by choosing an appropriate 
size of the step width we search for a relation for R such 
that one MCS corresponds to a real-time interval, in the 
sense of LD. 






le+OC 



100000 



10000 




1000 



0.1 
trial step width R 
FIG. 1. Characteristic time versus trial step width for a 
MC simulation. The solid line is fitted, yielding r ~ R^^. 



MC methods calculate trajectories in phase space fol- 
lowing a master equation which describes the coupling of 
a system to the heat bath. Hence, only the irreversible 
part of the dynamics of the system is considered ||l^ — 
there is no precession of the moments since no equation 
of motion is solved during the simulation. Nevertheless, 
in the following we will argue that the exact knowledge 
of the movement of the single moments is not necessary 
in order to describe the effects of thermal activation in 
an ensemble of systems under the following conditions: 
(i) the relevant time scales are larger than the precession 
time tp of the moments, (ii) we consider the high damping 
limit of the LLG equation where the energy dissipation 
during one cycle of the precession is considerably large 
so that the system relaxes (to the local energy minimum) 
on the same time scale tr ~ tr,. 
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FIG. 2. Configuration in phase space [Sx, Sy) of an ensem- 
ble of 20 particles following from a LD simulation for q = 1: 
see text for details. 

In Fig. we present the time evolution of our system 
in phase space, {Sx,Sy), following from a simulation of 
the LLG equation for high damping, a = 1. We use 
AE/ksT — 8.2, a rather low temperature so that the 
characteristic time r for the escape from the local en- 
ergy minimum is of the order of lO^^s (see also Fig. 13). 



The spin-precession time is tp = 9 x 10~"'^^s here. The 
simulation starts close to the local energy minimum with 
Sx — Sy — 0, Sz — 1 and the solid line shows the tra- 
jectory of one moment over a time interval of At — tp. 
The 20 points are the positions of an ensemble of 20 mo- 
ments after the same time. As one can see, the moments 
show no significant precession (the precession of an undis- 
turbed moment, i. e. without relaxation and fluctuations 
is indicated by the circle around the energy minimum at 
Sy = 0, Sx ^ 0.22). The small dots represent 1000 states 
of the ensemble for t < 6 x tp. Altogether, Fig. ^ demon- 
strates that in the high damping case already after time 
periods of only a few tp the moments are uncorrelated 
and the ensemble reaches a local equilibrium configura- 
tion (remember that the time scale to leave the local 
equilibrium is much larger here so that Fig.y shows only 
the local short-time equilibration, not the escape from 
the local energy minimum). 

We will show that this high-damping scenario can also 
be simulated by a MC simulation and we will now de- 
rive a relation for R in order to quantify the MC time 
step. The intention is to compare the fluctuations which 
are established in the MC technique within one MCS 
with the fluctuations within a given time scale associ- 
ated with the linearized LLG equation. Close to a local 
energy minimum one can write the energy, given that 
first order terms vanish as 



E K^ Eq + — 2_^ AijSiSj, 



(3) 



where the Si are the variables representing small devi- 
ations from equilibrium. In our system, for Bx = we 
find equilibrium along the z axis, leading to variables 
Sx and Sy. The energy increase AE associated with 
fluctuation in Sx and Sy is AE sa ^{AxxS^ + AyySy), 
with Axx ~ Ayy = 2dV -\- fisBz- Rewriting the LLG 
equation in the linearized form, Sx — LxxSx + LxySy, 
Sy = LyxSx + LyySy, It has bccu shown 
correlation function takes the form 



{S,{t)Sj{t'))^fi,,S,Jit~t'). 



that the 



(4) 



Dirac's 6 function is here an approximation for expo- 
nentially decaying correlations on time scales t — t' that 
are much larger than the time scale of the exponen- 
tial decay tr. The covarianz matrix fiij can be calcu- 
lated from the system matrices A and L as |1J] ^ij = 
— kBT{LikA^J + LjkA^^). For our problem a short cal- 
culation yields ^xx = fJ'yy = 2fcBT ^ "L — . Integrating 
the fluctuating magnetisation Sx (t) over a finite time in- 
terval At, Eq. Q takes the form 

(sl) = ^^xxAt = 2fcBr "^ A^, (5) 

(1 + a^)/is 

representing the fluctuations of Sx averaged over a time 
interval At. 



Next, we calculate the fluctuations (S*^) during one 
MCS of a MC simulation. This is possible if we as- 
sume that all magnetic moments are initially in their 
equilibrium position. For our MC algorithm described 
above the probability distribution for trial steps with 

step width r = JS^ + SI is p, = iy/W^ / {2-k R^) . 

The acceptance probability within a heat bath algorithm 
isp,(r) = l/(l-hexp(AS(r2)/fcBT)), where AEir"^) can 
be taken from Eq. |[ Hence, for the fluctuations within 
one MC step it is: 



{SD = 27r 



'rdryp,(r)p.(r) = :^+0(i?4) (6) 



where the last line is an expansion for small R. By equal- 
izing the fluctuations within corresponding time intervals 
we find the relation 



R' 



20kBTaj 



At. 



(7) 



Note, from our derivation above it follows that one time 
step At must be larger than the intrinsic time scale t^ of 
the relaxation. This means - as already mentioned above 
- that the Monte Carlo method can only work on time 
scales that are much larger than any microscopic time 
scale of a precession or relaxation (to local equilibrium) 
of the moment. 

In principle, equation m gives the possibility to choose 
the trial step for a MC simulation in such a way that 
IMCS corresponds to a real time interval, say At = 
10^ ^■^s. However, there are of course restrictions for pos- 
sible values of R, like R < 1. Also, R should not be too 
small since then a Monte Carlo algorithm is inefficient. 
Therefore, either one has to choose such a value for At 
so that R takes on reasonable values (these will usually 
be of the order of 10~^^s) or one uses a reasonable con- 
stant value for R, say 0.1, and uses Eq. to calculate 
At as the real time interval corresponding to IMCS. In 
the following we use the first method since it turns out 
to be very efficient to change R with temperature. How- 
ever, we confirmed that the other method yields the same 
results. 
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FIG. 3. T versus inverse temperature: comparison of IHD 
asymptote, LD simulation, and MC simulation. 



To test the validity of our considerations we performed 
MC simulations with an algorithm using a trial step ac- 
cording to Eq. with Ai « 6 x 10~^^s (the inverse value 
of 7, in other words the time in the LLG equation is 
rescaled by 7). For Fig. we set a = 1 and compare 
the data for t{T) following from our MC simulation with 
results from LD simulations and with the intermediate 
to high damping (IHD) asymptote MM, namely 



2ttujo 



,/3(Vb-V2) 



2TrujQ 



.AE/keT 



(8) 



where uq and Qq are the saddle and damped saddle an- 
gular frequencies which have been defined in Eqs. (21) 
and (22) of Ref. |6) explicitly. W2 is the well angular fre- 
quency for the deeper of the two potential wells and is 
defined in Eq. (20) of Ref. Q. All have been defined 
in terms of the coefficients of the truncated Taylor se- 
ries representation of the energy equation described in 
detail in section V of Ref. B, (particularly Eqs. 59-64). 
For the purpose of comparison with MC and LD simula- 
tions, we consider one escape path only, e^^^°~^^\ where 
/3 = V/kBT and Vq — V2 is the energy described by Eq. 
(62) of Ref. Q. For our purposes, /3(Vo — V2) may be 
represented by AE/kBT. The validity condition for the 
IHD formula is aAE/ksT > 1 where AE/ksT > 1 
which have been satisfied in all cases represented here. 

From Fig. y it is clear that the LD data agree with the 
asymptote above. For higher temperatures the asymp- 
tote is no longer appropriate. Here, the numerical data 
for T tend to zero for T — + cxd as one expects. The MC 
data deviate slightly and are roughly 10% larger. How- 
ever, considering the fact that to the best of our knowl- 
edge this is the first comparison of a " real-time MC sim- 
ulation" with LD simulations and asymptotic formulae, 
the agreement is remarkable - especially taking into ac- 
count the simple form of Eq.^ underlying our algorithm 
and also that there is no adjusted parameter in all our 
calculations and formulae. 
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FIG. 4. r versus damping constant: comparison of IHD 
asymptote, LD simulation, and MC simulation. 

Since we expect that our MC procedure leads to a 
high damping limit we also tested the a-dependence of r. 
Fig. H shows the corresponding data for the same param- 
eter values as before and AE/T = 3.3. The figure shows 
that the MC data converge to the IHD formula and to 



the data from LD simulation for large a. Even the small 
10% deviation of the MC data mentioned before (Fig. |^) 
vanishes in the limit of larger a. 

To summarize, we discussed the conditions under 
which a comparison of LD with a MC process appears 
to be possible. Considering a simple system of isolated 
single-domain particles we derived an equation for the 
trial step width of the MC process so that one step of 
the MC algorithm can be related to a certain time inter- 
val. Testing this algorithm we found excellent agreement 
with data from LD simulation as well as with interme- 
diate to high damping asymptotes for the characteristic 
times of the magnetisation reversal. Even, though our al- 
gorithm was derived only for the special system which we 
consider here, we belive that the arguments we brought 
forward might the fundament even for the MC simulation 
of more complicated systems, especially systems consist- 
ing of interacting magnetic moments. 
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